The sleep EEG envelope is a novel, neuronal firing-based human biomarker

Sleep EEG reflects voltage differences relative to a reference, while its spectrum reflects its composition of various frequencies. In contrast, the envelope of the sleep EEG reflects the instantaneous amplitude of oscillations, while its spectrum reflects the rhythmicity of the occurrence of these oscillations. The sleep EEG spectrum is known to relate to demographic, psychological and clinical characteristics, but the envelope spectrum has been rarely studied. In study 1, we demonstrate in human invasive data from cortex-penetrating microelectrodes and subdural grids that the sleep EEG envelope spectrum reflects neuronal firing. In study 2, we demonstrate that the scalp EEG envelope spectrum is stable within individuals. A multivariate learning algorithm could predict age (r = 0.6) and sex (r = 0.5) from the EEG envelope spectrum. With age, oscillations shifted from a 4–5 s rhythm to faster rhythms. Our results demonstrate that the sleep envelope spectrum is a promising biomarker of demographic and disease-related phenotypes.

The sleep EEG is a continuous signal reflecting ongoing electrical activity in the brain, and its spectrum reflects the relative contribution of different frequencies to the final waveform. In contrast, the envelope of the sleep EEG estimates the instantaneous amplitude of the signal (typically after filtering to frequencies of interest), and its spectrum estimates the periodicity of all band-limited activity. In other words, the envelope spectrum estimates the typical rhythm at which signal amplitude at certain frequencies waxes and wanes (Fig. 1).
In contrast there has been no systematic research describing the periodicity, generating mechanism or reallife correlates of the EEG envelope, in spite of the fact that as a mathematical function of the oscillations from which the ordinary sleep EEG spectrum is calculated it has the theoretical potential to be an equally promising biomarker with potential incremental validity. On June 1, 2021 we searched PubMed and ScienceDirect with the search terms "eeg envelope spectrum", "eeg envelope psd" and "eeg envelope power". Our search returned no relevant papers. We also screened the first 100 Google Scholar hits with these search terms, but also found no relevant papers. Nevertheless, we are aware of some previous studies which do not explicitly assess the spectrum of the EEG envelope but which are still of interest to this field. An early paper 42 described that sleep spindles followed each other in periods slightly exceeding four seconds, corresponding to a hypothetical 0.25 Hz envelope oscillation. Two studies 43,44 calculated PSD from short windows, smoothed the resulting power estimates and relied on the spectral analysis of the resulting signal to Figure 1. The principle of EEG envelope spectrum analysis. (A) Shows a simulated EEG signal, consisting of the sum of a 2 Hz sinusoid modulated by a 0.2 Hz carrier frequency, a 12 Hz sinusoid modulated by a 1 Hz carrier frequency, and pink noise. Overlain blue and red lines show the instantaneous amplitude or envelope (modulus of the Hilbert transform) of the delta (1-4 Hz) and sigma (10)(11)(12)(13)(14)(15)(16) Hz) frequency ranges, respectively. (B) Shows the power spectral density of the original signal (left) and the delta (middle) and sigma (right) envelopes. Note that the carrier frequencies are accurately recovered from spectral analysis of the envelopes (with some impurities due to added noise and the fact that the modulus of the Hilbert transform of a modulated signal is not fully sinusoidal). The spectrum of the envelope reveals periodic fluctuations in the amplitude of higher-frequency activities.

Results
The envelope reflects cortical neuronal firing. In Study 1, we used invasive human EEG data from epileptic patients. We correlated multiple-unit activity (MUA) measured by cortex-penetrating microelectrodes with the envelope of EEG signals measured on the adjacent cortical surface. We found that in all frequency bands and across the entire cortical mantle, the envelope of the surface signal reflected firing patterns, with a typical average magnitude-squared coherence value of 0.15-0.2. The pattern of coupling was different as a function of frequency range (Fig. 2).
In the delta through the theta range, MUA was lowest during or slightly after the peak of the envelope, also reflected by the fact that the largest cross-correlations were negative and observed when a zero to slightly negative MUA delay was added. In the sigma through gamma ranges, MUA was highest during the ascending phases of the envelope, also reflected by the fact that the largest cross-correlations were positive in case of a positive MUA delay. The alpha range exhibited an intermediate pattern with only a few IME channels reaching significance.
Findings in individual patients are available in the Supplementary Data. See also Fig. 3B for an illustration of envelope-MUA coupling. We note that in 3 patients the envelope-MUA coupling was absent or restricted to very specific channels. As data issues (problems with synchronization in case of an absent coupling, and poor MUA data quality in case of both absent and spatially restricted coupling) is a possible explanation for this pattern, we re-analyzed coupling excluding these three patients. Results were virtually identical even in this case ( Supplementary Fig. S1).
The envelope spectrum is stable within individuals. For Study 2, we used scalp EEG data from healthy volunteers to calculate the spectrum of the envelope using a novel colliding window method (see "Methods" and Fig. 3A). This method ensured optimal data availability in case of artifacts (Fig. 3C), which is a particular concern for the estimation of very low frequency spectral power which requires long windows. We found characteristic spectral peaks at ~ 0.05 Hz (20 s periodicity) for all frequency bands, ~ 0.25 Hz (4 s periodicity) for several frequency bands, in particular the sigma band corresponding to sleep spindles, and ~ 1 Hz mainly for the beta and gamma bands (Fig. 3D).
The ordinary sleep EEG spectrum is known to have a trait-like quality by being stable within individuals, but varying between individuals 2,4 . In the absence of multiple recordings from participants, in Study 2 we assessed the trait-like nature of the envelope spectrum by calculating even-odd reliabilities (intraclass correlation coefficients between the spectral densities calculated separately from the even or odd numbered sampling windows of the same individual) and split-half reliabilities. Split-half reliability is a method frequently used in psychometric literature. It refers to the similarity of the first and second halves of a non-instantaneous measurement (most typically, the scores derived from the first and second halves of a psychological questionnaire) 54 . In our current interpretation, split-half reliabilities were defined as the Pearson correlations between the spectral densities calculated separately from only the first and last 50% sampling windows of the same individual.
Based on this analysis, the envelope spectrum was highly trait-like in NREM sleep and moderately so in REM sleep. The mean reliability of the envelope EEG spectrum (pooled across channels, frequency bands and envelope frequencies) was 0.886 (even-odd, SD = 0.085) and 0.819 (split-half, SD = 0.14) for NREM and 0.513 (even-odd, SD = 0.157) and 0.684 (split-half, SD = 0.141) for REM. No clear trend was seen by envelope frequency (Fig. 4). The reliability of the mid-frequency EEG bands (alpha and low sigma) was the highest, falling off towards both low and high frequencies ( Supplementary Fig. S2), while no clear trend was seen by scalp channel (Supplementary Fig. S3).
The envelope spectrum reflects age and sex, but not general cognitive ability. The envelope spectrum of the sleep EEG was significantly associated with demographic variables, but not with general cognitive ability (Figs. 5,6). This association was the strongest between the NREM envelope PSD and age. Older age was generally associated with a loss of low-frequency oscillations in the power of NREM EEG frequencies, often with an increase of higher-frequency oscillations. Specifically, a reduced oscillation of low delta power at ~ 0.5 and ~ 1.6-1.8 Hz, but increased oscillation at 1-1.5 Hz; a reduced ~ 0.25 Hz oscillation of theta, alpha, sigma and beta power with an increased 0.5-1 Hz oscillation of theta and sigma power was seen. In REM sleep, a general tendency for increased low-and high-frequency power oscillations and a corresponding decrease at ~ 0.5-1 Hz was seen, but this only reached statistical significance in the high delta, alpha and low sigma frequency bands.
Male sex was associated with a lower amplitude of ~ 0.05-0.1 and ~ 0.5-1.5 Hz NREM low sigma power oscillations, but a higher amplitude of power oscillations of the same frequency band at ~ 0.25-0.5 Hz and > 1.75 Hz. Male sex was also associated with a lower amplitude of < 0.75 Hz and > 1. 75 Hz, but a higher amplitude of ~ 1 Hz beta power oscillation irrespective of sleep state.
General cognitive ability was not significantly associated with the envelope spectrum of either NREM or REM sleep EEG.
Multivariate models. The relationship between human phenotypes and single biological markers, such as single genetic polymorphisms or individual features of brain morphology is usually modest. However, multivariate models using a large number of such biological markers as independent variables are able to capture the additive, independent contribution of each single marker to reach a much more substantial correlation between the totality of biological markers and phenotypes 56 . Therefore, beyond demonstrating correlations between single spectral features of the sleep EEG envelope and age, sex and intelligence, we set out to investigate the relationship between these features using multivariate models.  www.nature.com/scientificreports/ tion sample (N = 22), we only used the first 100 PSD bins as these exhibited the largest bivariate correlations with phenotypes, and we ran models separately by channel. The total predictive validity of the envelope spectrum towards each phenotype was expressed as the correlation between predicted and actual values in the validation sample (predictive accuracy) (Fig. 7). Using NREM sleep, age could be predicted with reasonable accuracy (r mean = 0.616, r SD = 0.151, the prediction accuracy for sex was lower but still substantial (r mean = 0.447, r SD = 0.138), but the correlation between predicted and actual IQ was low (r mean = 0.151, r SD = 0.178). Using the REM sleep envelope, age could be predicted with moderate accuracy (r mean = 0.502, r SD = 0.156), but this was not the case for sex (r mean = − 0.019, r SD = 0.224) and IQ (r mean = 0.092, r SD = 0.112). (All means and SDs are across channels.) In case of IQ, elastic net models frequently failed to converge due to the low correlation between PSD values and this phenotype.
Relationship to respiratory rhythms. Low-frequency fluctuations in the EEG signal could theoretically be affected or contaminated by the respiratory cycle, which also occurs with sub-second periodicity. In order to investigate to what extent this occurs, we used the recordings of 29 participants containing a pair of respiration channels to estimate to what extent respiratory activity is correlated with the course of EEG envelopes. We estimated the (1) magnitude-squared coherence between each EEG band envelope and respiratory activity (2) the modulation index between respiratory activity and EEG band envelopes. We performed these calculations with 100 s windows of artifact-free data with 50% overlaps between windows and compared statistics to those calculated from 1000 random surrogates to estimate statistical significance in each participant. Like in other analyses, we transformed the resulting empirical p-values to z-values, averaged them across participants and transformed them back to p-values before application of the Benjamini-Hochberg correction of false discovery rate.
In line with a previous study 43 we found no coupling between EEG envelopes and respiratory activity. Neither coherence between respiratory activity and EEG envelopes nor their modulation index was ever significantly higher than in surrogates, illustrated by an almost perfectly circular phase histograms of EEG envelope amplitudes as a function of respiration phase (Supplementary Figs. S4, S5). Thus, our findings confirm that EEG amplitude fluctuations occur largely independently from low-frequency respiratory rhythms, and thus the EEG envelope is not a respiratory artifact.
Relationship to the ordinary EEG spectrum. In contrast to previous studies investigating the spectrum of the EEG signal, we investigated the spectrum of the envelope of the signal. Although there is a mathematical relationship between these two measures, it is unclear how well they are correlated at the individual level. In order to investigate this, we calculated correlation coefficients between the individual mean values of the envelope spectrum (calculated as described above, in all frequency bands and both vigilance states) and the spectrum of the log-transformed absolute signal PSD (calculated with the Welch method, using 4 s epochs with an 50% overlap, in both vigilance states, yielding a spectral resolution of 0.25 Hz). The envelope spectrum was considered in the 0.01-1 Hz range and the signal spectrum was considered in the 0-48 Hz range.In both vigilance states, the correlation between the signal PSD and the envelope PSD was generally low, suggesting that the two measures are largely statistically independent. A full illustration of all correlations are provided on Supplementary Figs. S6 and S7.
However, some correlations were present: • NREM low delta envelope: positive correlation with < 2 Hz signal PSD at 0.01-0.2 Hz, negative correlation at > 0.5 Hz. • NREM theta and alpha envelope: positive correlation with ~ 8-12 Hz (both band envelopes) and > 15 Hz (theta only) signal PSD at 0.2-0.3 Hz. • REM theta, alpha and low sigma envelopes (the latter probably an extension of alpha): positive correlation with mid-frequency (~ 5-25 Hz, with an even wider range for low sigma) signal PSD at < 0.1 Hz, a negative correlation with the same signal PSD frequencies at 0.2-0.8 Hz. • REM beta envelope: a positive correlation with 1-5 Hz signal PSD at 0.3-0.6 Hz.
Correlations between the NREM low delta envelope and the low delta PSD likely reflect that in participants with more slow wave activity, this activity tends to fluctuate in rather long cycles. The other observations need replication and specific experimental designs to enable a physiological interpretation, but overall they highlight that fluctuations in EEG power-especially in the mid-frequencies-captured by the envelope spectrum are associated with individual differences in signal PSD beyond the frequency range from which the envelope was calculated. This effect was especially pronounced in REM sleep.

Discussion
In our study, we aimed to describe the sleep EEG envelope in detail and compare its characteristics to the ordinary sleep EEG spectrum to assess its viability as a biomarker. Overall, our study demonstrates that the sleep EEG envelope shares many of the properties of the ordinary sleep EEG: it reflects neuronal population firing, it has characteristic oscillation frequencies, it is highly individually stable and varies between individuals; and it is associated with demographic characteristics.
It has been shown in previous human invasive EEG studies that sleep oscillations recorded either from the cortex or from the scalp closely reflect the rhythmic ensemble firing of neuron populations. For instance, slow waves 11,12 , sleep spindles 57 and the wakeful alpha rhythm 15 as field potentials are all associated with waxing and waning patterns in local neuronal firing. We observed a similar pattern for the envelope as well. MUA was www.nature.com/scientificreports/ significantly suppressed when low-frequency activity was high: specifically, the lowest MUA was observed during the maximum of low delta activity and slightly after the maximum of high delta and theta activity. Curiously, the opposite pattern (increased MUA during periods of reduced low-frequency oscillations) was less typical. This phenomenon may reflect the rhythmic suppression of neuronal firing during slow oscillations 11,12,58 , which contain ensembles of low frequencies up to the alpha range 59 . Although such neuronal down-states are generally followed by up-states containing high frequency rhythms 52,60 , the fact that down-states are generally more prominent 12,58 may specifically result in a specific association between the presence of low-frequency activity in the ECoG and reductions in neuronal firing in nearby cortex. For high frequencies (low sigma through gamma, with alpha being an intermediate range), an opposite pattern was seen: MUA was maximal when oscillations in Illustrates the log-transformed envelope spectra. All data was z-transformed by frequency band to eliminate mean differences. The frequency axis is shown on a log scale to enhance the low frequency ranges which are of particular interest. Note spectral peaks at ~ 0.05-0.06 Hz, ~ 0.25 Hz and ~ 1 Hz, the latter most prominent in the beta range. www.nature.com/scientificreports/ these frequencies were gaining in power, possibly reflecting the role of cortical neuronal assemblies in recruiting these oscillations. Next, we used scalp EEGs for healthy volunteers to establish further properties of the sleep EEG envelope spectrum. We found that, similarly to the ordinary spectrum, the envelope spectrum was also characterized by higher powers at lower frequencies. In line with previous reports on slow wave and sleep spindle periodicity 42,43 we found two characteristic peaks: one at ~ 0.05 Hz (20 s period, most prominent for slow rhythms), and another at ~ 0.25 Hz (4 s period, most prominent for faster rhythms). These frequency peaks were less prominent in REM sleep than in NREM.
Previous reports have established that the sleep EEG spectrum is fingerprint-like with a high intra-individual stability [1][2][3][4][5] , which is the result of genetic regulating factors [6][7][8][9] . Although our ability to fully replicate this finding in the envelope spectrum was limited by the absence of multiple recordings and genetically informative data, we could establish that when comparing spectra from the same individual across the two halves of the night or across even and odd numbered sampling windows, reliability was very high for NREM (> 0.8) and reasonably high for REM (> 0.5), with remarkably similar reliability values across all but the lowest frequencies.
The reliability of the sleep EEG envelope spectrum renders it a potential marker of stable individual differences, such as demographic variables, psychological traits or pathological conditions. In a quantitative test of this hypothesis, we found that higher age was associated with reductions in the 0.25 Hz rhythmicity of high delta through beta rhythms. A relative increase in the ~ 1 Hz rhythm of sigma-frequency oscillations, an additional increase in the very low frequency rhythms of low sigma and beta oscillations, as well as a relative reduction of low-frequency and a relative increase of high-frequency low delta oscillations was also seen. These resultstogether with findings from invasive EEG recordings-can be interpreted as a systematic loss of the medium-scale temporal organization of rhythmic neuronal firing as a function of ageing. Notably, envelope spectra calculated from NREM were much more associated with age than REM spectra, highlighting the functional importance of this vigilance state for ageing-related phenomena 19,20,41,61 .
Sex was associated with a single EEG envelope feature: low-frequency rhythmicity of the NREM beta rhythm was reduced in males, while high-frequency rhythmicity was higher. The significance of the ~ 1 Hz rhythm suggests that beta rhythms show stronger coupling to slow waves in males, however, the functional importance of this finding is currently unknown.
Although intelligence was found to be associated with multiple sleep EEG spectral features 25,62 , we found no evidence that it is also associated with the rhythmicity of sleep EEG oscillations.
We used a learning algorithm to perform multivariate predictions of age, sex and intelligence based on the sleep EEG envelope spectra. As expected based on the reliability of spectra, much better predictions could be made based on NREM than REM spectra. Age could be predicted with reasonable accuracy from NREM sleep envelope spectra (r ~ 0.6), although much more accurate predictors were previously constructed based on overall features of the sleep EEG 19 or the shape of NREM slow waves 20 Sex could be predicted from the NREM envelope spectrum with lower but still substantial accuracy (r ~ 0.45), although the predictive power of the REM   www.nature.com/scientificreports/ spectrum was much lower. Sex prediction based on the envelope spectrum underperforms relative to other predictors based e.g. on brain imaging [63][64][65] . However, we did not expect the envelope spectrum to be a particularly sexually dimorphic characteristic. The non-significant zero-order correlations between the envelope spectrum and intelligence could not be improved with the use of elastic net regression: models failed to converge on most electrodes and even with this method we found no association between the envelope spectral and intelligence. What biological process do amplitude fluctuations in the sleep EEG reflect? In Table 1 we provide a nonexhaustive list of known biological oscillations with periods at most on the minute scale. From this list, we had data about two prominent oscillations: the cardiac and the respiratory rhythm. Both oscillations could theoretically drive low-frequency EEG rhythms either through physiological mechanisms (for example, because neuronal firing depends on the availability of oxygenated blood and this is reflected in EEG rhythms) or through electrical artifacts detected by the EEG. However, based on non-significant magnitude-squared coherence and phase-amplitude coupling the respiratory rhythm appears to play a role in low-frequency EEG amplitude oscillations, and the cardiac rhythm is too fast to strongly influence all but the fastest envelope rhythms. Because our recordings did not contain data about other oscillating biological processes, we can only speculate about their role. With their characteristic 20-s periods, gastric rhythms 66 oscillate at a frequency strongly overlapping with characteristic envelope frequencies, rendering EEG envelope oscillations a promising potential marker in the study of brain-viscera interactions. Other known biological rhythms are not strong candidates to be the driver of or to be coupled with EEG amplitude oscillations due to differences in their characteristic frequencies. In sum, the precise biological mechanism creating periodic fluctuations in the amplitude of EEG rhythms remains unknown and its discovery is a major task of further studies into this phenomenon.
In sum, our study revealed that the periodicity of amplitude fluctuations in the sleep EEG, reflected by the envelope, is a promising human biomarker. In an invasive study, it was found to be associated with fluctuations in neuronal firing. In a study of healthy volunteers, it was found to be a highly reliable individual marker, somewhat sexually dimorphic and especially strongly associated with ageing. While we showed that envelope fluctuations reflect fluctuations in neuronal firing, why these fluctuations take place (and why they change with ageing) requires further study.
Our work has a number of limitations. First, using a single IME per patient we were only able to record neuronal firing from a very limited cortical area. EEG recorded on the adjacent cortical surface is likely the summation of neuronal activity in a more extended area, consequently, the correlation between MUA and the envelope was not particularly strong. Second, we had only a single night of measurement from healthy individuals, resulting in a within-night, rather than a more optimal across-night estimation of envelope reliability.

Methods
Participants. Study 1. Sleep electrophysiological data from 13 patients undergoing presurgical electrophysiological evaluation for drug-resistant epilepsy were used. All interventions were approved by the Hungarian Medical Scientific Council and the ethical committee of the National Institute of Clinical Neuroscience. Clinical procedures were not biased for scientific purposes. All patients gave informed consent in line with the Declaration of Helsinki.

Study 2.
We used data from 176 healthy participants (mean age 29.8 years, SD 10.66 years, range 17-69 years; 95 males) from a multi-center database of the Max Planck Institute of Psychiatry (Munich, Germany) and the Psychophysiology and Chronobiology Research Group of Semmelweis University (Budapest, Hungary) 20,68 was used in this retrospective study. We used participants with available cognitive test scores (Raven's Advanced Progressive Matrices, the Culture Fair Test and/or the Zahlenverbindungstest [a trail making test]). Test scores were always expressed as IQ scores with a population mean of 100 and a standard deviation of 15, and if multiple tests were available from a single participant, the scores were averaged (see the first publication of the dataset 68 for details). Study procedures were approved by the ethical boards of Semmelweis University, the Medical Faculty of the Ludwig Maximilian University or the Budapest University of Technology and Economics. All participants were volunteers who gave informed consent in line with the Declaration of Helsinki. According to semi-structured interviews with experienced psychiatrists or psychologists, all participants were healthy, had no history of neurologic or psychiatric disease, and were free of any current drug effects, excluding contraceptives in females. Consumption of small habitual doses of caffeine (maximum two cups of coffee until noon), but no alcohol, was allowed. Six male and two female participants were light-to-moderate smokers (self-reported), and the rest of the participants were non-smokers. Further details about participant selection criteria and study protocols can be found in the studies reference above.
Electrophysiology. Study 1. Patients underwent electrophysiological recordings using implanted laminar microelectrodes (IME) and subdural grid and strip electrodes, from which only grids were analyzed (ECoG). Detailed descriptions of these methods are described elsewhere 11,57,69 . In brief, IMEs contain 24 serially referenced contacts on a cortex-penetrating pin spaced evenly at 150 µm, capable of detecting extremely local intracortical electrical activity, including neuron population firing, which is represented by high-frequency data (300-5000 Hz) from this source. Multiple-unit activity (MUA), an index of local neuronal population firing, was calculated by rectifying raw data and filtering it with a 20 Hz low-pass filter, according to standard procedure 11,57 . ECoG was recorded with a sampling frequency/precision of either 2000 Hz/16 bit or 1024 Hz/16 bit depending on the individual patient, and always with a contralateral mastoid reference.
We manually selected seizure-free data with adequate signal quality (indicated by the absence of continuous, broad-frequency artifacts) from all patients. Sleep staging for the selected ECoG data was performed visually on a 20 s basis based on standard criteria 70  www.nature.com/scientificreports/ Figure 6. Correlations between the REM envelope spectrum and age, sex and general cognitive ability (IQ). Colored lines represent correlation coefficients by scalp channel. Color codes indicate scalp region, with individual channels from the same region shown with the same color. Black horizontal lines show the threshold of conventional (p = 0.05) significance. Colored dots (with color coding identical to lines) above the lines indicate a statistically significant correlation after FDR correction on the corresponding channel. www.nature.com/scientificreports/ EEG channels with a full polysomnography setup (including EOG and EMG), we restricted our scoring to the identification of NREM sleep (regardless of stage) and the separation of it from other sleep states and wakefulness, based on the presence of slow waves and spindles. REM sleep, which is difficult to detect using our setup, was not analyzed in Study 1. Artifacts were excluded from ECoG data on a 4 s basis using visual inspection. Only artifact-free data from NREM sleep was considered for further analysis. For analysis, we selected the ECoG channel closest to the IME without epileptiform activity. For the IME, we treated data from poor-quality channels (based on visual inspection) as missing data.

Study 2.
All participants underwent all-night polysomnography recordings for two consecutive nights, and data from the second night was used for all analyses. Scalp EEG electrodes were applied according to the 10-20 system and referenced to the mathematically linked earlobes. Impedances were kept at < 8kΩ. EEG was sampled at 250 Hz for 115 participants, 249 Hz for 29 participants and 1024 Hz for 15 participants, always resampled at 250 Hz. Sleep EEG was visually scored on a 20 s basis according to standard criteria 70 . A visual scoring of artifacts was also performed on a 4 s basis. EEG preprocessing was implemented in Fercio's EEG (©Ferenc Gombos, Budapest, Hungary). Further details about the technical details of the sample can be found in the first publication of this dataset 68 . For analyses, we pooled all NREM epochs instead of analyzing N2 and SWS separately. This was motivated by observations that sleep depth within NREM is continuous rather than categorical 71 and concerns that age-related and sleep depth-related changes in slow wave activity may be confounded in our demographically heterogeneous sample 60 .
ECoG-MUA coupling. For each segment and for each frequency band, we estimated the coupling between ECoG envelope and the MUA by calculating (1) the normalized cross-correlation of the two signals implemented with the xcorr() MATLAB function, allowing lags in the [− 1 1] second range; (2) the magnitudesquared coherence between the two signals at 0.1 Hz intervals between 0.1 Hz and 1 Hz, implemented with the mscohere() MATLAB function; and (3) the coupling of the amplitude of ECoG envelopes to MUA phases. For this last analysis, we first used the phase angle of the Hilbert transform to estimate the instantaneous phase of the MUA signal. Next, we z-transformed the ECoG envelope signal along the time dimension to standardize amplitude across segments. Finally, we calculated the mean standardized ECoG envelope amplitude (expressed in within-segment SD units) concomitant to MUA data in each of 12 equally spaced phase bins of 30 degrees each. For each patient, we averaged each of the three statistics across all segments to generate a mean value. We used this method to estimate phase-amplitude coupling because traditional methods 72 , only estimate the preferred phase and overall significance of coupling, whereas we aimed to calculate a more fine-grained estimate. We note, however, that our method is theoretically closest to the Modulation Index 73 , except we estimate the statistical significance of each histogram bin individually instead of relying on a single, Shannon entropy-based estimate of omnibus significance.
Statistical significance calculation. We estimated the statistical significance of coupling statistics by comparing results to surrogates obtained from random EEG segments. For this, we matched each 20-s ECoG envelope segment with a randomly selected artifact-free NREM MUA segment, calculated cross-correlation, coherence and phase-amplitude coupling and finally an average value across all segments. We performed this analysis 1000 times to generate a null distribution of coupling statistics. An empirical p-value was assigned to each statistic based on actual data, defined as the proportion of surrogate-based statistics more distant from zero.
We calculated unweighted means of all comparable statistics across patients. Similarly, we transformed p-values into standard normal deviates (z-scores) and averaged them across patients, similarly to Fisher's method of averaging logarithmized p-values 74 This approach is more conservative and different from ordinary meta-analysis in that it doesn't add weights to patients based on the amount of data available and it doesn't increase power over what was originally available in individual patients, so effects which fall short of significance in individual patients do not become significant when data is pooled. Effectively, the alternative hypothesis of this method is that coupling is significantly different from zero in each patient, while in a standard meta-analysis it would be that coupling is significantly different from zero when data from all patients is pooled.
Finally, average standard normal deviates were transformed back to p-values and subjected to correction for false discovery rate using the Benjamini-Hochberg method 75  www.nature.com/scientificreports/ band (cross-correlation), across frequencies by frequency band and IME channel (coherence), and across phase bins by frequency band and IME channel (phase-amplitude coupling).
Colliding window method and spectrum smoothing. Because fluctuations in the envelope of the EEG signal are expected to take place on a much longer timescale than fluctuations in the signal itself, very low frequencies of the envelope spectrum are of particular interest, but their estimation is only possible with sampling windows much longer than those used to estimate the ordinary power spectrum. This introduces a particular problem when dealing with artifacts. In case of the ordinary power-spectrum, which is estimated using many sampling windows each only a few seconds long, the loss of a few sampling windows due to the presence of artifacts only results in the loss of a comparatively small fraction of the total signal. In case of the envelope spectrum, however, totally discarding a 100-s sampling window due to a presence of a relatively short artifact may result in an unacceptable amount of signal loss. Therefore we used a colliding window method (Fig. 2, Panel A) to deal with artifacts. When the 100-s windows sampling the signal in 20 s steps encountered a segment marked as an artifact, they were progressively shortened to end before the artifact, until a minimum sampling window length of 20 s was reached. At this point, the sampling window skipped the artifact segment and re-started at its original 100-s duration afterwards. PSD from the shortened windows was calculated and used as usual, but PSD estimates of the frequencies below 1/L Hz were discarded and in the calculation of the average PSD data from this window was under-weighted by 1*L/100 (L in both cases refers to the length of the window in seconds). In order to avoid over-sampling of data before artifacts, all envelope signals were sampled both in the forward and backward direction, starting the 100 s windows from the beginning and the end of recordings, respectively. The colliding method ensured a minimum signal loss of 20 s instead of 100 s in case of artifacts. The resulting average envelope PSDs were smoothed using the Savitzky-Golay method with a 10-degree polynomial, 10-base log-transformed to normalize variances for linear statistics and z-transformed across frequencies within participant and channel to eliminate the effect of between-participant differences in raw EEG signal voltage. Participants with abnormal PSDs (based on visual inspection) on any channel in any frequency band were removed from analyses concerning that frequency band (N = 1-3 participants per frequency band).
Envelope frequencies up to 2 Hz (that is, fluctuations in EEG amplitude with up to two cycles per second) were considered for analysis. Figure 2 illustrates the colliding window process, the amount and temporal position within the night of available artifact-free data and the average spectra. Detailed individual spectra are available in the Supplementary Data. Envelope reliability. Even-odd reliability was computed by calculating the average PSD for each individual twice, using even and odd numbered sampling windows separately. Because sampling windows were up to 100 s long and overlapped by 20 s steps, only every fifth sampling window was used to avoid non-independent data. The reliability of the PSD in each frequency band, on each channel and at each frequency was estimated by the intraclass correlation coefficient (implemented as Pearson's correlation coefficient with pooled standard deviations) between the two measurements. For split-half reliability, we also calculated the average PSD for each individual twice, using the first and last 50% of all available sampling windows separately. Because the intraclass correlation coefficient is sensitive to mean differences and we expected mean signal voltage to systematically change between the first and second halves of the night, we computed split-half reliability using the ordinary Pearson correlation instead. Although reliability is generally defined as the square root of the correlation between repeated measurements because they are both expected to be equally affected by unreliability 76 , we used the more conservative and more easily interpretable unsquared coefficients.
Multivariate analysis. For multivariate predictions, we used elastic net regression implemented in the MAT-LAB lasso() function. Elastic net regression is an iterative learning algorithm which seeks to maximize the predictive value of a large number of potentially correlated predictors by introducing a penalty term for complexity. Elastic net regression is able to fit reliable models in samples where OLS regression would be underdetermined given the large number of predictors and the small sample size. Technical descriptions 77,78 and practical implementations 79,80 , including in sleep EEG analysis 20 are available in the literature. We used fivefold crossvalidation and an L1-L2 regularization mixture set at alpha = 0.5 for elastic net regression models. All envelope spectral values between 0.01-1 Hz from all spectral bands (800 variables in total) were used as predictors and age, sex (here treated as a continuous variable 81 ) and IQ were used as dependent variables. These models were fitted independently using data from each electrode (18*3 = 54 models in total). One eighth (N = 22) of the sample was retained as a validation sample, and the models were trained on the remaining participants (N = 154, including the cross-validation samples used to ensure robust regression coefficients). The models resulting from training were used in the fully independent validation sample to check performance. The validation sample was selected by ordering participants by the values of the dependent variable and taking every 8th individual to ensure maximal variance.

Data availability
Supplementary data is available on Zenodo at https:// doi. org/ 10. 5281/ zenodo. 55953 41. Due to limitations described in the ethical permit of this study (especially pertaining to patient data), raw EEG data is available upon reasonable request to the corresponding author. All original code has been deposited at Zenodo at https:// doi. org/ 10. 5281/ zenodo. 55953 41 Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.